Adaptive Thouless- Anderson-Palmer approach to inverse Ising problems with 

quenched random fields 



- 1—1 

X 



c3 



PACS numbers: 02.50.Tt, 02.30.Zz, 75.10.Nr 



I. INTRODUCTION 



Haiping Huang and Yoshiyuki Kabashima 
Department of Computational Intelligence and Systems Science, 
Tokyo Institute of Technology, Yokohama 226-8502, Japan 
(Dated: May 7, 2013) 

The adaptive Thouless- Anderson-Palmer equation is derived for inverse Ising problems in the 
presence of quenched random fields. We test the proposed scheme on Sherrington-Kirkpatrick, 
Hopfield, and random orthogonal models and find that the adaptive Thouless-Anderson-Palmer 
approach allows surprisingly accurate inference of quenched random fields whose distribution can 
be either Gaussian or bimodal, compared with other existing mean-field methods. 

in ■ 

The inverse Ising problem has been intensively studied in statistical physics and computational biology in the past 
few years [H-Q. Such studies are of huge practical and theoretical relevance. On one hand, the advent of techniques 
for multi-electrode recording and microarray measurement produces high-throughput biological data Unveiling 
the biological mechanism underlying these experimental data poses a challenging computational problem. In the 
inverse Ising problem, one tries to construct a statistical mechanics description of the original system directly from 
the data, and it provides a promising tool for dimensional reduction in modeling vast amounts of biological data [(| ■ 
On the other hand, for guaranteeing the reliability of the obtained description, it is also necessary to examine the 
reconstruction performance of the inverse algorithms numerically and/or analytically by utilizing artificial data that 
I ' are generated from a variety of known Ising spin models [U, 0413] ■ 

In general, the experimental data are described by M independent samples {er 1 , er 2 , . . . , cr M } in which er is an 
TV-dimensional vector with binary components (<7; = ±1) and N is the system size. The Ising model provides the 
O | least structured model to match the statistics of the experimental data as 



cn: Pisi - (CT) = z(ib) exp 

o ■ 



(ij) 1 



(1) 



where (ij) denotes all distinct spin pairs and the partition function Z(h, J) depends on iV-dimensional fields and 

^vq ' — —-dimensional couplings. These fields and couplings are chosen to yield the same first and second moments 
(magnetizations and pairwise correlations, respectively) as those obtained from the experimental data. The inverse 
f-^ temperature (3 = l/T is absorbed into the strength of fields and couplings. 

. Based on magnetizations and correlations, the inference of fields and couplings of the Ising model is a computa- 
— ' tionally hard problem especially for large systems. However, one can resort to mean-field methods, such as naive 
mean-field (nMF) [11], Thouless-Anderson-Palmer (TAP) equation 0, Sessak-Monasson (SM) expansion [l2|, and 
Bcthe approximation (BA) [l3|-[l5j , to get an approximate solution to the inverse problem with computationally fea- 
sible costs. Previous investigations have mostly focused on the inference of the coupling vector, whereas the inference 
5h ■ error of fields has been less analyzed. In fact, external fields represent intrinsically preferred directions of {a - ;}, which 
are also very important for understanding information processing in real neuronal networks [l], [T(| and gene interaction 
networks 01 and for predicting protein structures from sequence data @, [l7[- Therefore, an accurate estimation of 
external fields is also highly desirable. 

To this end, we propose the adaptive Thouless-Anderson-Palmer (adaTAP) approach for the inverse Ising problem 
and establish the framework on the basis of Gibbs free energy and Gaussian approximation. Surprisingly, adaTAP 
yields a very accurate estimation of external fields, although some other mean-field methods are more competitive 
in predicting couplings. We confirm the efficiency of adaTAP on three kinds of mean-field models: the Sherrington- 
Kirkpatrick (SK) model [3, the Hopfield [H model and the random orthogonal model (ROM) [20(; other existing 
mean-field inverse algorithms are also compared. 

The outline of this paper is as follows. The adaptive TAP approach to the inverse Ising problem with quenched 
random fields is derived in Sec. [TTJ In Sec. Mil extensive numerical simulations are carried out to test the inference 
performance of adaTAP on Hopfield model, SK model and ROM. The comparison with other existing mean-field 
methods is also made and discussed. Concluding remarks are given in Sec. [IV] 
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II. ADAPTIVE TAP APPROACH 



For the Ising model defined in Eq. (p}, we write the magnetization-dependent free energy (also termed Gibbs free 
energy) as 



G(m) = -h T m + Extr e j^m - ln^] . 



(2) 



where the Lagrange multiplier vector 8 is introduced to fix magnetizations at all sites to their thermal expectation 
values, i.e., rrii — (er,:). h T denotes the transpose of a vector h. The notation Extr stands for the extremum with respect 
to the corresponding parameters (6 here). The exact evaluation of the partition function in Eq. is computationally 
difficult for a large system. However, one can resort to mean-field approximations. We adopt the following strategy. 
First, each coupling is multiplied by a real number I € [0, 1], and the Gibbs free energy can then be expressed by 



G(m) = G(m,l = 1) = f 
Jo 



8G(m,l) 



= D= / dl y —^+G(m,l = 0) 

o 91 

G g (m, / = 1) - G g (m, I = 0) + G(m, I = 0), 



G(m, I) = h T m + Extr e j 9 T m - In ^ 
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G g (m, I) = h 1 m + Extr e W 1 m - -tr(AC 



In / d(Te -^ T (A-u) CT +e r <, 



(3a) 
(3b) 

(3c) 



where Cy = {(JiUj), and we have used Gaussian statistics for the binary spins with expectation constraints, i.e., 
(<7iCTj) g = (cjCTj) Is j ng which are enforced by a symmetric matrix A. Here, tr(A) denotes the trace of a matrix A. 
For simplicity, we assume A is a diagonal matrix A = diag(Ai, . . . , A^r) whose diagonal terms are determined via the 
extremization of the corresponding Gibbs free energy. The Gaussian approximation makes the computation of the 
partition function tractable. This scheme is also called the expectation consistence approximation (2lj and was applied 
to derive the message-passing algorithm for the perceptron learning problem (22|. Conventional Plefka expansion 23 1 
truncates the power series expansion of G(m, I) to second order in /, but Eq. (|3a| contains terms of all orders. Note 
that the third term in Eq. (J3a|) is the Gibbs free energy of non-interacting Ising spins at fixed magnetizations and can 
be easily evaluated. The final expression for the Gibbs free energy reads, 



G(m) ~ -im T Jm - h T m + ^2'H{m l ) + i lndct(A - J) 

i 



(4) 



where l-L(mi) 



l+m, i l+m, 
2 111 2 



1 2 m ' In 1 2 mi and A follows the extremization condition of Eq. ([3c|l with 1 = 1, 



(A - J),: 1 = 1 - m\ 



(5) 



Equilibrium values of magnetizations are determined by m oq 
A quick calculation gives the self-consistent equation for m, 



argmin m G(m) and the free energy F — min m G(m) 



tanh 



m j 



A, 



1 



(G) 



which is exactly the adaptive TAP equation first introduced in Refs. [24], [23 for the Ising model. Eq. © can also 
be derived under other mean-field approximations (26l - [28j . The third term inside the square bracket of Eq. ([6]) forms 
the Onsager correction term which requires no prior knowledge of the coupling statistics, playing an important role 
in inferring external fields. Aj in Eq. © is a function of {mj determined by Eq. ([S]). The fixed point of the self- 
consistent equation gives m oq . We remark here that Eq. ^ can be reduced to the normal TAP equation obtained 
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from a high-temperature expansion of the Gibbs free energy 23, 29, 30], i.e., the third term inside the square bracket 
of Eq. becomes — Xy (1 — rn j)^ij (Onsager reaction term) through high-temperature expansion. 

To obtain the inference equations for couplings, we use the identity HC = I 0, IU HH where H is the Hessian 
matrix of the Gibbs free energy ify = g^.g m and C is the connected correlation matrix whose entries are Cy = 

(o~iO~j) — rriimj. I is an identity matrix. Magnetizations and correlations are already given by the experimental data 
in the inverse Ising problem. Finally, the inference equation reads 

J« = -(C- 1 ) ij + m i (B-% (7) 

for i ^ j, where By = [Xij] 2 which expresses how large the change of mi is given a small perturbation to Aj and 
we define x = (A — J) -1 - The expression for Bij is derived by using the Sherman-Morrison formula [33]. A small 
perturbation AAj to Aj will lead to a corresponding change of m, according to Eq. (O, which is described by the 
following equation: 

AA-y 2 - 

[X- 1 + AA^ 1 = Xii 7 J = 1 - K + Am,) 2 , (8) 

1 -t- L^l\jXxj 

where the Sherman-Morrison formula is used in the first equality and the notation AAj means that only j-th diagonal 
term of matrix AA is non-zero and equal to AAj. Noting that both Ami and AAj are small, one can obtain 
= = 2^7- [Xij] 2 by using Eq. (0 once again. After couplings are reconstructed, external fields are inferred as 

hi = tantrum;) - ^2 Jio m j + m i [Aj - - — — ^ J . (9) 

To predict the coupling, we need to solve the adaTAP equation Eq. ([5|). An iterative scheme is proposed as follows. 
Step 1. Let t = and initialize Jy = — (C _1 )y, A,; = (1 — mf ) _1 for all (ij) and i respectively. 

Step 2. At t, set t' = 0, A*' = ° = A*, x = (A* - J*) -1 . 

Step 2.1. t' <- t' + 1, update A*' = A*'^ 1 + A A, where A A; = jz^js ~ ^ for au After update of each Af , x 

y old AAy old 

needs to be updated simultaneously as Xfc? w — Xki ~ i+aa x°* d ^ er ^ ve( i by using the Sherman-Morrison 
formula. 

Step 2.2 Until |AA.;| < ea for all i, then assign A* — A , x l = X and go to Step 3. Otherwise, if t' < <J nax , 
go to Step 2.1, else return UN-CONVERGED. 

Step 3. t «-t+l, update 4 = -(C-%- + m^B^y where By = ^- [x*~ 1 f, until |Jf - J^r 1 ] < tJ for all (ij), 
then go to Step 4. Otherwise, if t < t max , go to Step 2, else return UN-CONVERGED. 

Step 4. Infer according to Eq. © for all i. 

In step 2.1, the step size AA^ for updating Aj can be derived by using Eq. (|S|) and the Shermon-Morrison formula, 
which gives xu ~ i+aa*^ = ^ ~ ^ n ^ e iterative scheme, we set the parameters i max = t' max = 1000, and 
£A = e./ = 10~ 4 . In practice, we find that both A and J in our simulations shown below converge in tens of steps 
when the temperature is not very low. The computational complexity of this iterative scheme is dominated by the 
inverse of the matrix (e.g., C or B), keeping the same order as that of other mean-field methods. 

To compare performances of different mean-field inverse algorithms, we define the inference error for couplings and 
fields, respectively, as 



f VfJ* - J 1 

N(N-l)^ Kv 13 



true \ 2 



1/2 

(10a) 



±-Y( h *-htn 2 



N 

i 

where J*- (h*) is the inferred coupling (field) and J*™° (/i' ruo ) is the true one 



1 1/2 

(10b) 
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III. NUMERICAL SIMULATIONS 



We evaluate the inference performance of the adaTAP approach on three mean-field models with either Gaussian 
distributed or bimodal distributed random fields. For the SK model, each entry of the coupling matrix are indepen- 
dently drawn at random from a Gaussian distribution with zero mean and variance 1/N. In the Hopfield model, the 
coupling is constructed according to Hebb's rule, i.e., Jy = jj Ylu=i wnere P random Gaussian patterns 
are stored in the network, are independent Gaussian random variables with zero mean and unit variance. We also 
test our method on ROM whose coupling matrix is constructed as J = O t ~DO where O is an orthogonal matrix chosen 
with the Haar measure [20, |34j- D = diag(Ai, . . . , Ajy) and A follows a distribution p(A) = aS(X — 1) + (1 — a)S(X + 1). 

To collect the data of magnetizations and correlations, we use the exact enumeration on small-size systems of 
N = 15, which produces noise-free data for predicting the underlying parameters. In this case, M — 2 N . Inference 
results of adaTAP on these three tested models are compared with those obtained by other mean-field methods. 
For comparison, we briefly describe the other four existing mean-field methods for the inverse Ising problem. The 
couplings between spin i and j (i ^= j) are inferred as follows: 



J. 



J, 



nMF 

ij 

TAP 



2(C" 1 ) 



J™ = 



-1 - y/1 ~ 8mim,j(C 



7" 



MF 



Tind 



L. t Lj — Cfj 



jB A = -tanh- 1 



2(C 



rT — (aij - bij) 



(11a) 
(lib) 

(11c) 
(lid) 



where J, 



ind 



(l+C, J ) 2 -(m i +m 3 ) 
(l-C (j ) 2 -(m,-mj)' 



, (l « „,.,. • Oij = ^l + AL t L 3 {C-%, bij = y /(a ij -2m i m j (C-i) ij ) 2 -4(C-% and 

Li = 1 — mf. After couplings are inferred, fields can be predicted using the following equations: 

^ MF =tanh- 1 (m,)-^ J$ 



qMF. 



7 TAP„ 



hJ AP = tanlT 1 (rrn ) - J ij • • 
+ miJ2(^ AP ni-m 2 j ), 



hf A = tanh 1 (m^) — tanh 1 (tijf(mj,mi,tij)) 



(12a) 

(12b) 
(12c) 



where = tanh Jjj and f(x, y, t) 



_ l-t 2 -y/(l-t 2 ) 2 -4t(x-yt)(y-xt) 



2t(y-xt) 



Since SM expansion has large inference errors for 
predicting fields even when considering up to the third order in the small correlation expansion [l2j . we would not 
show its field inference performances for the temperature range we consider. 

We first examine the inference performance of adaTAP on mean-field models, where quenched random fields are 
drawn independently at random from a Gaussian distribution with zero mean and variance a\. As displayed in 
Fig. Q] (a) for the Hopfield model, adaTAP shows slightly better performance than the TAP approach in coupling 
constructions, whereas the SM expansion has the best performance at high temperatures and the BA has the best 
one at low temperatures. Regarding field inference, adaTAP performs much better than other methods in the entire 
temperature ran ge u nder consideration. However, if we incorporate an effective self-coupling (diagonal weight) Ja — 
1 _ L m i — (C^ 1 )a [llj into the inference equation (|12a[) . nMF with diagonal weights (nMFdw) will achieve the nearly 
same accuracy with adaTAP in predicting external fields, although adaTAP still gives a bit lower inference error. This 
also holds for the other two mean-field models. Note that nMF without diagonal weights definitely gives a highest 
inference error among all mean-field methods compared here. As the temperature becomes sufficiently low, adaTAP 
ceases to converge within i max or t' ma:K , thus becoming unable to predict couplings and fields. To infer a model with 
quenched random fields, TAP and BA will also have no solution at low enough temperatures. We also performed 
simulations with a larger <j\ (e.g., a\ = 0.1), and it is observed that the field inference performance deteriorates 
and adaTAP fails to converge at a higher temperature for some samples compared to the case with a smaller field 
variance. Fig.[T](b) shows inference results for ROM with the random orthogonal coupling matrix. Although adaTAP 
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FIG. 1: (Color online) Inference performances of adaTAP on Hop field, random orthogonal and SK models with Gaussian 
distributed random fields, compared with those obtained by other existing mean-field methods. Magnetizations and correlations 
used to infer fields and couplings are calculated through exact exhaustive enumeration on networks of size N = 15. Each data 
marker is the average over 20 random realizations for which g\ — 0.01. (a) Results for the Hopfield model with P = 3. (b) 
Results for ROM with a = 0.6. (c) Results for the SK model. 

behaves slightly worse than TAP for inferring couplings, it produces surprisingly accurate estimates of external fields 
in the entire temperature range in Fig. [T](b). Note that the inference accuracy obtained by other mean- field methods 
(except nMFdw) can be further improved by at least one order of magnitude by using adaTAP when the random 
fields are Gaussian distributed. For coupling inferences of the SK model (see Fig.[T](c)), the performance of adaTAP 
lies between those of nMF and TAP, while the SM expansion gives a more accurate prediction than other methods at 
high temperatures. 

Fortunately, the superiority of adaTAP for field inference is also true when the random field is bimodal distributed, 
i.e., Ph(h) = p5{h — ho) + (1 —p)8(h + ho). Its performance is shown in Fig. [5] with p — 0.6, ho = 0.3. The improvement 
of the field prediction by adaTAP is evident in this case, even compared to nMFdw. In adaTAP, we still assume 
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FIG. 2: (Color online) Inference performances of adaTAP on Hopfield, random orthogonal and SK models with bimodal 
distributed random fields, compared with those obtained by other existing mean-field methods. Magnetizations and correlations 
used to infer fields and couplings are calculated through exact exhaustive enumeration on networks of size N = 15. Each data 
marker is the average over 20 random realizations for which ho = 0.3, p = 0.6. (a) Results for the Hopfield model with P = 3. 
(b) Results for ROM with a = 0.6. (c) Results for the SK model. 



zero diagonal couplings, however, the third term inside the square bracket of Eq. ([6]) provides an adaptive Onsager 
correction to the nMF approximation, playing the same key role with diagonal couplings in inferring external fields. 
In this adaptive manner, lower inference error of fields and couplings is achieved compared to nMFdw. Interestingly, 
adaTAP can even perform better than TAP in predicting couplings for certain ranges of temperatures in this case. 

All the models investigated so far are of the fully connected type. For examining the capability to deal with another 
extreme of sparsely connected networks, the proposed scheme is also tested on a tree model. Our proposed scheme 
performs better than other methods except BA which is exact on a tree and gives very accurate inference both on 
the couplings and fields. The result is shown in Fig. [3] A tree of size = 22 is constructed, such that each node 
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FIG. 3: (Color online) Inference performances of adaTAP on a tree (N = 22), compared with nMFdw and BA (exact method 
on a tree model). Each data marker is the average over 10 random realizations. 

inside the tree has degree equal to 3, and to mimic an infinite Bethe lattice, we generate the external fields for the 
boundary spins as hi — hi + J2kedi\j hk->i, where j is the only spin inside the tree connected to the boundary spin i 
and the cavity field /ifc_>j is randomly chosen from a population dynamics for an infinite Bethe lattice (3oT | . Couplings 
and fields (hi for the boundary spins) for the tree follow Gaussian distributions Af(0, 1) and A/70, 0.01) respectively. 
Magnetizations and correlations are calculated by using susceptibility propagation algorithms [16[ . However, in real 
applications, for example, a typical neuronal network of size around ./V = 100 is not strongly diluted with an exact tree 
structure, therefore, our method is expected to still give good estimates of external fields. To confirm this point, we 
test adaTAP on a diluted SK model, where each non-zero Gaussian distributed coupling is present with a predefined 
probability pd- The Gaussian distribution has zero mean and variance 1/c with c = pdN. As shown in Fig. 21 the 
adaTAP still performs better than other mean-field methods (including nMFdw) in field inference, which is much 
more evident when random fields are bimodal distributed. 

Although adaTAP, TAP, and the BA will have no solution in inferring mean-field models with quenched ran- 
dom fields at low temperatures, adaTAP does outperform other existing mean-field methods compared here to infer 
quenched random fields if it converges, as shown for a wide range of temperatures on the Hopfield, random orthogonal 
and (diluted) SK models. We conclude that the power of adaTAP for inverse Ising problems resides in its remarkable 
accuracy in predicting external fields, especially for the case where there is a single dominant state in phase space. 

IV. CONCLUSION 

In summary, we propose the adaTAP approach for inverse Ising problems and show its striking performance for 
inferring external fields in mean-field models. As far as the field inference is concerned, adaTAP is rather satisfactory, 
compared to other mean-field methods. Furthermore, an accurate inference of external fields in the Ising model is 
able to provide us with insights into the mechanism underlying high-throughput data either coming from biological 
experiments or from large databases [l], [U, [TtJ ■ The proposed adaTAP approach for the inverse Ising problem is 
expected to have applications in real data analyses (e.g., neural data, or sequences in the protein databases), in 
combination with other mean-field methods. 
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